CTB Modeling
Rita Ludwig, John Flournoy
2019-02-13 12:22:27
We will fit a multi-level non-linear (on the parameters) regression model to the data from the sample. The regression takes the form:
\[ \text{sooner_choice}=\frac{\text{endowment_later}\cdot(b^{t_0}d^k\text{pratio})^{\frac{1}{a-1}}}{1+\text{pratio}\cdot(b^{t_0}d^k\text{pratio})^{\frac{1}{a-1}}} \] where \(\text{endowment_later}=20\).
In R, this looks like the following formula for the population-level parameters (i.e., the fixed effects):
sooner_choice ~ ( endowment_later * (b^t0 * d^k * pratio)^(1/(a-1)) ) / ( 1+pratio * (b^t0 * d^k * pratio)^(1/(a-1)) )We start by importing the data and ensuring that it contains all of the information for each trial we need. For example, we need to add columns that specify the delay period for each trial.
mctb_data_ <- utils::read.table(file.path(data_dir_base, 'mCTB_kit/data.txt'),
header = TRUE,
na.strings = '.')
mctb_data_l <- stats::reshape(mctb_data_,
varying = paste0('c', 1:24),
v.names = 'c',
timevar = 'budget_number',
times = 1:24,
direction = 'long')
mctb_design <- read.table(file.path(data_dir_base, 'mCTB_kit/instrument_details.txt'),
header = TRUE)
mctb_data <- merge(mctb_data_l, mctb_design,
by = 'budget_number',
all.x = TRUE, all.y = TRUE)
head(mctb_data)## budget_number subject_id c id endowment_soon endowment_later
## 1 1 1 6 1 19 20
## 2 1 2 6 2 19 20
## 3 1 3 1 3 19 20
## 4 1 4 1 4 19 20
## 5 1 5 1 5 19 20
## 6 1 6 3 6 19 20
## sooner_date_weeks delay_weeks
## 1 0 5
## 2 0 5
## 3 0 5
## 4 0 5
## 5 0 5
## 6 0 5We will also take this opportunity to examine missing reponses.
n_missing_responses <- dplyr::summarize(group_by(mctb_data, subject_id),
n_missing = sum(is.na(c)),
p_missing = n_missing/n())
plot(n_missing_responses$subject_id, n_missing_responses$p_missing,
type = 'h', xlab = 'Participant ID', ylab = 'Proportion missing responses',
ylim = c(0,1))
knitr::kable(filter(n_missing_responses, p_missing > 0),
digits = 2,
caption = paste0('Proportion missing by participant. Total N with missing data = ',
sum(n_missing_responses$p_missing > 0)))| subject_id | n_missing | p_missing |
|---|---|---|
| 110 | 3 | 0.12 |
| 217 | 1 | 0.04 |
| 270 | 20 | 0.83 |
| 395 | 24 | 1.00 |
| 613 | 20 | 0.83 |
| 624 | 20 | 0.83 |
| 655 | 6 | 0.25 |
| 1007 | 2 | 0.08 |
| 1068 | 2 | 0.08 |
Here, we define some trial-level predictors. t0 defines whether the soonest choice is today or also in the future. We modify k so that it is in terms of days, rather than weeks. We also calculate the price ratio, which is just the ratio of the later reward and sooner reward.
mctb_data$t0 <- as.numeric(mctb_data$sooner_date_weeks == 0)
mctb_data$k <- 7 * mctb_data$delay_weeks
mctb_data$k_w <- mctb_data$delay_weeks
mctb_data$pratio <- mctb_data$endowment_later / mctb_data$endowment_soonHere we calculate the possible value choices the participant saw for each trial.
for(j in 1:6){
mctb_data[paste0('soon_', j)] <- (20/as.vector(mctb_data$pratio)) - (j-1) * (20/as.vector(mctb_data$pratio)) / 5
mctb_data[paste0('late_', j)] <- (j-1) * 4
}
mctb_data$soon_6 <- 0Finally, based on the particpant’s response, 1-6, we create a column that contains the amount of money they choose to receive at the soonest possible time. If they choose option 6, to receive the full $20 at the later time, they always receive $0 at the sooner time.
#Create a new variable with a value that corresponds to the response option
#chose in `mctb_data$c`
mctb_data$sooner_choice <- apply(
mctb_data, 1, #for every row in mctb_data
function(a_row){
#write the name of the column to get the choice value from
choice_col <- paste0('soon_', a_row['c'])
#get the choice value
sc <- a_row[choice_col]
#if it's null, return NA, otherwise return the choice value
return( ifelse(is.null(sc), NA, sc))
})Exclusion criteria are specified:
If participants miss greater than or equal to 50% of items on a personality scale, they will be coded as missing a score for that scale and not be included in analyses involving that scale. Due to the nature of the CTB task, participants who miss more than one item per timeframe pair will be coded as missing and excluded from analysis. Finally, participants who did not self-report household income will be excluded from analysis.
As we saw above, some participants are missing more than one response on the CTB task. No one did not report income. Some participants failed to answer >50% of items on a personality scale.
scale_data <- read_csv(file.path(data_dir_base, 'JRPdataOSF.csv'))## Parsed with column specification:
## cols(
## .default = col_double(),
## annualincome = col_character()
## )## See spec(...) for full column specifications.scale_data_missing <- select(scale_data, contains('score'), subID) %>%
filter_at(.vars = vars(-subID), any_vars(is.na(.)))
missing_ids <- full_join(filter(n_missing_responses, n_missing > 0),
scale_data_missing,
by = c('subject_id' = 'subID'))
mctb_data <- anti_join(mctb_data, missing_ids,
by = 'subject_id')Total number of excluded IDs is 21, or 1.9%.
We use fixed-effects non-linear least-squares as a first pass on the data and to set starting values for the individual-level analyses that will assist the multi-level estimation.
sd_of_soonerchoice <- sd(mctb_data$sooner_choice, na.rm = TRUE)
#scale the outcome
mctb_data$sooner_choice_scaled <- mctb_data$sooner_choice / sd_of_soonerchoice
mctb_data$endowment_later_scaled <- mctb_data$endowment_later / sd_of_soonerchoice
a_start <- 0.90
b_start <- 1
d_start <- 0.999
nls_ctb_model <- nls(sooner_choice ~
( endowment_later * (b^t0 * d^k * pratio)^(1/(a-1)) ) /
( 1+pratio * (b^t0 * d^k * pratio)^(1/(a-1)) ),
mctb_data,
start = list(a = a_start,
b = b_start,
d = d_start))
knitr::kable(broom::tidy(nls_ctb_model), digits = 4)| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| a | 0.7451 | 0.0049 | 152.8705 | 0 |
| b | 0.9608 | 0.0045 | 213.8533 | 0 |
| d | 0.9986 | 0.0000 | 21015.3706 | 0 |
##STORE MODEL PARAMS##
alpha_nls <- summary(nls_ctb_model)$coefficients[1]
beta_nls<- summary(nls_ctb_model)$coefficients[2]
delta_nls <- summary(nls_ctb_model)$coefficients[3]
mctb_data$alpha_nls <- alpha_nls
mctb_data$beta_nls <- beta_nls
mctb_data$delta_nls <- delta_nlsWe can interpret these values as follows: \(a\) is considerably less than 1, indicating that generally people do not prefer the all-sooner or all-later options, but, rather, prefer to take some money sooner and some money later. The parameter \(b < 1\) indicates people discount even a bit more by a factor of \(b = 0.9608\) if they can get the sooner reward today, rather than waiting five weeks. In other words, for a constant delay of 5 weeks, they discount at a rate \(d=0.9986\) exponentiated by day then multiplied by \(0.9608\) if they get the sooner reward today (the sooner discount for a five week delay is \(0.9986^{5\times7} = 0.9523\), but if the sooner reward is today, it’s \(0.9608 \times 0.9986^{5\times7} = 0.915\)).
If we apply this to the full equation above, we can find the sooner payout that, combined with the portion of the payout received later, is equivalent to a $20 later payout for a participant with these parameter values (if the sooner payout is received today, for a price ratio of, say, 20/16 = 1.25):
round((sooner_choice_calc <- (20 * (beta_nls*delta_nls^(7*5)*1.25)^(1/(alpha_nls-1)) )/
( 1 + 1.25 * (beta_nls*delta_nls^(7*5)*1.25)^(1/(alpha_nls-1)) )),2)
round((amount_received_later <- 20-sooner_choice_calc*20/16), 2)
round(total_amount_received <- sooner_choice_calc + amount_received_later, 2)## [1] 6.79
## [1] 11.51
## [1] 18.3\[ 6.79 = \frac{20\cdot(0.9608^{1}\cdot 0.9986^{7\times5}\cdot1.25)^{\frac{1}{0.7451-1}}}{1+1.25\cdot(0.9608^{1}\cdot 0.9986^{7\times5}\cdot1.25)^{\frac{1}{0.7451-1}}} \]
So, these population parameters suggest that someone will be just as happy getting $6.79 today and $11.51 (for a total of $18.3) later as they would to receive the full $20 later.
Another way of expressing the discounting variable is in terms of a monthly discount rate:
(a_monthly_rate <- delta_nls^(-30) - 1)## [1] 0.04275648Individual-level parameters
Finding the cases that have no meaningful variability in their responses. We exclude trials where you could receive $20 sooner or $20 later.
odd_cases <- group_by(mctb_data, subject_id) %>%
filter(soon_1 != late_6) %>%
summarize(
no_interior = case_when(
any(c %in% 2:5) ~ FALSE,
all(c %in% c(1,6)) ~ TRUE,
TRUE ~ NA
),
no_variance = case_when(
all(c == 1) ~ TRUE,
all(c == 6) ~ TRUE,
TRUE ~ FALSE
))
print(paste0('Do all cases that have no variance also have no interior solutions? ',
all(odd_cases$no_interior[odd_cases$no_variance])))## [1] "Do all cases that have no variance also have no interior solutions? TRUE"#View(odd_cases)
table(odd_cases[,2:3])## no_variance
## no_interior FALSE TRUE
## FALSE 533 0
## TRUE 201 345This indicates that 345 participants always answer either sooner or later exclusively (excluding the $20 now versus $20 later trials). 533 participants use the interior options, and 201 participants just switch between taking everything sooner, or everything later.
mctb_nlslist <- nlsList(
sooner_choice_scaled ~
( endowment_later_scaled * (b^t0 * d^k_w * pratio)^(1/(a-1)) ) /
( 1+pratio * (b^t0 * d^k_w * pratio)^(1/(a-1)) ) | subject_id,
data = select(mctb_data,
t0, k_w, pratio, subject_id,
sooner_choice_scaled, endowment_later_scaled),
start = c(a=alpha_nls,b=beta_nls, d=delta_nls),
control = list(maxiter = 5000, minFactor = 1/2^30, tol = 1e-1))## Warning: 492 errors caught in numericDeriv(form[[3L]], names(ind), env). The error messages and their frequencies are
##
## step factor 4.65661e-10 reduced below 'minFactor' of 9.31323e-10
## 9
## number of iterations exceeded maximum of 5000
## 14
## singular gradient
## 100
## Missing value or an infinity produced when evaluating the model
## 369converged <- lapply(mctb_nlslist,
function(amod) !is.null(amod))
converged_df <- data_frame(subject_id = as.numeric(names(converged)),
nlslist_converged = unlist(converged),
nlslist = 1)
mctb_data_diagnostics <- left_join(
left_join(mctb_data, odd_cases,
by = 'subject_id'),
converged_df[,-3],
by = 'subject_id')Take a look at responses from those with convergent versus non-convergent nls models. The first three panels are convergent, and the last three are non-convergnet.
ggplot2::theme_set(ggplot2::theme_minimal())
plot_mctb_responses(
filter(mctb_data, subject_id %in%
c(converged_df$subject_id[converged_df$nlslist_converged][1:3],
converged_df$subject_id[!converged_df$nlslist_converged][1:3])),
adj.pratio = TRUE)
Use the nlslist object to start up the linear mixed effects model:
mctb_nlme.nlslist <- nlme.nlsList(model = mctb_nlslist,
random = b + d ~ 1,
groups = ~subject_id,
method = 'ML')
summary(mctb_nlme.nlslist)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: sooner_choice_scaled ~ (endowment_later_scaled * (b^t0 * d^k_w * pratio)^(1/(a - 1)))/(1 + pratio * (b^t0 * d^k_w * pratio)^(1/(a - 1)))
## Data: select(mctb_data, t0, k_w, pratio, subject_id, sooner_choice_scaled, endowment_later_scaled)
## AIC BIC logLik
## 51574.3 51631.44 -25780.15
##
## Random effects:
## Formula: list(b ~ 1, d ~ 1)
## Level: subject_id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## b 0.13088065 b
## d 0.02301261 0.704
## Residual 0.57988744
##
## Fixed effects: list(a ~ 1, b ~ 1, d ~ 1)
## Value Std.Error DF t-value p-value
## a 0.9454517 0.000794510 24815 1189.9811 0
## b 0.9190270 0.004295194 24815 213.9663 0
## d 0.9829883 0.000715799 24815 1373.2737 0
## Correlation:
## a b
## b -0.027
## d -0.016 0.621
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.578681e+00 -2.167673e-01 -1.581553e-05 2.787358e-01 4.000113e+00
##
## Number of Observations: 25896
## Number of Groups: 1079qplot(sample = residuals(mctb_nlme.nlslist)/sd(residuals(mctb_nlme.nlslist)),
geom = c('qq', 'abline'),
main = "Q-Q plot for conditional residuals")
qplot(sample = ranef(mctb_nlme.nlslist)$b/sd(ranef(mctb_nlme.nlslist)$b),
geom = c('qq', 'abline'),
main = "Q-Q plot for random b")
qplot(sample = ranef(mctb_nlme.nlslist)$d/sd(ranef(mctb_nlme.nlslist)$d),
geom = c('qq', 'abline'),
main = "Q-Q plot for random d")
qplot(x = predict(mctb_nlme.nlslist), y = residuals(mctb_nlme.nlslist),
geom = c('jitter', 'smooth'), width = .05, height = .05)## Warning: Ignoring unknown parameters: height## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
qplot(x = predict(mctb_nlme.nlslist), y = sqrt(abs(scale(residuals(mctb_nlme.nlslist)))),
geom = c('jitter', 'smooth'), width = .05, height = .05)## Warning: Ignoring unknown parameters: height## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
plot(mctb_nlme.nlslist, subject_id ~ resid(.))
plot(mctb_nlme.nlslist, sooner_choice_scaled ~ fitted(.) | subject_id)mctb_nlme.nlslist_coef_df <- as.data.frame(coef(mctb_nlme.nlslist))
mctb_nlme.nlslist_coef_df$subject_id <- as.numeric(rownames(mctb_nlme.nlslist_coef_df))
mctb_data_model_info <- left_join(mctb_data_diagnostics,
mctb_nlme.nlslist_coef_df,
by = 'subject_id')
mctb_model_results <- distinct(mctb_data_model_info,
subject_id, a, b, d,
no_interior, no_variance, nlslist_converged)
write_csv(mctb_model_results,
path = file.path(data_dir_base,
'JRP_CTB_model_paramsOSF.csv'))mctb_pars_and_flags <- distinct(mctb_data_model_info,
a, b, d,
no_interior, no_variance, nlslist_converged,
subject_id)
flag_summaries <- group_by(mctb_pars_and_flags,
no_interior, no_variance) %>%
summarize(a = mean(a, na.rm = T),
b = mean(b, na.rm = T),
d = mean(d, na.rm = T),
n = n())
knitr::kable(bind_rows(flag_summaries,
as.data.frame(t(fixef(mctb_nlme.nlslist)))), digits = 4)| no_interior | no_variance | a | b | d | n |
|---|---|---|---|---|---|
| FALSE | FALSE | 0.9455 | 0.9211 | 0.9838 | 533 |
| TRUE | FALSE | 0.9455 | 0.9282 | 0.9851 | 201 |
| TRUE | TRUE | 0.9455 | 0.9105 | 0.9805 | 345 |
| NA | NA | 0.9455 | 0.9190 | 0.9830 | NA |
for(avar in list(list('b', .05), list('d', .005))) {
print(ggplot(mctb_pars_and_flags,
aes_string(x = avar[[1]],
group = 'nlslist_converged',
fill = 'nlslist_converged')) +
geom_vline(xintercept = 1) +
geom_histogram(position = position_dodge(),
aes(y = ..count..),
binwidth = avar[[2]]) +
facet_grid(no_interior ~ no_variance,
labeller = labeller(
no_interior = as_labeller(c('TRUE' = 'No interior choices',
'FALSE' = 'Has interior choices')),
no_variance = as_labeller(c('TRUE' = 'No choice variance',
'FALSE' = 'Has choice variance'))),
as.table = TRUE) +
labs(x = paste0('Parameter ', avar[[1]]),
fill = 'Single-subject\nNLS model\nconverged?') +
theme_minimal())
}
Participants who have very low d values:
plot_small_d_sids <- unlist(distinct(filter(mctb_data_model_info,
d < .94, !no_variance),
subject_id))
for(i in unique((1:23 -1) %/% 6)){
last_index <- ifelse(i*6 < length(plot_small_d_sids),
(i+1)*6,
length(plot_small_d_sids))
first_index <- i*6 + 1
print(plot_mctb_responses(filter(
mctb_data_model_info,
subject_id %in% plot_small_d_sids[first_index:last_index]))
+theme_minimal())
}
Parameter-behavior correlation check
behavior_summary <- group_by(mctb_data, t0, subject_id) %>%
summarize(sooner_choice = mean(sooner_choice)) %>%
gather(key, value, -subject_id, -t0) %>%
unite(key, key, t0) %>%
spread(key, value) %>%
mutate(mean_choice = (sooner_choice_0 + sooner_choice_1)/2,
diff_choice = sooner_choice_1 - sooner_choice_0)
compare_df <- left_join(mctb_model_results, behavior_summary)## Joining, by = "subject_id"plot(compare_df$mean_choice, compare_df$d)
plot(compare_df$diff_choice, compare_df$b)
The above plots indicate a reasonable mapping of individual-level parameter values to the observed behavior.
Plotting all the choices
theme_set(theme_minimal())
mctb_sample_ids <- distinct(mctb_data, subject_id) %>%
mutate(plot_group = (1:n()-1) %/% 5)
for(g in 0:max(mctb_sample_ids$plot_group)){
mctb_data_sample <- left_join(
filter(mctb_sample_ids, plot_group == g),
mctb_data, by = 'subject_id')
aplot <- plot_mctb_responses(mctb_data_sample)
print(aplot)
}
Check sensitivity to letting \(\alpha\) vary over individuals
#this model shows poor convergence
mctb_nlme.nlslist_a <- nlme.nlsList(model = mctb_nlslist,
random = a + b + d ~ 1,
groups = ~subject_id,
method = 'ML',
control = nlmeControl(msVerbose = TRUE,
niterEM = 100,
pnlsTol = 1.1),
verbose = 2)## 0: 119821.65: 1.72729 1.86088 3.52593 -1.24060 -12.2423 -25.7977
## 1: 119821.65: 1.72729 1.86088 3.52593 -1.24060 -12.2423 -25.7977
## 2: 119821.65: 1.72729 1.86088 3.52593 -1.24060 -12.2423 -25.7977
## 3: 119821.65: 1.72729 1.86088 3.52593 -1.24060 -12.2423 -25.7977
##
## **Iteration 1
## LME step: Loglik: -24990.93, nlminb iterations: 3
## reStruct parameters:
## subject_id1 subject_id2 subject_id3 subject_id4 subject_id5 subject_id6
## 1.727291 1.860884 3.525932 -1.240604 -12.242275 -25.797727
## Beginning PNLS step: .. completed fit_nlme() step.
## PNLS step: RSS = 5412.874
## fixed effects: 0.8425702 0.9113531 0.9825489
## iterations: 5
## Convergence crit. (must all become <= tolerance = 1e-05):
## fixed reStruct
## 0.03503821 12.64544046
## 0: 115425.25: 0.380387 1.20965 2.81685 -0.914147 -0.897170 -12.8678
## 1: 115425.25: 0.380223 1.20944 2.81665 -0.914152 -0.897172 -12.8678
## 2: 115425.25: 0.380223 1.20944 2.81665 -0.914152 -0.897172 -12.8678
## 3: 115425.25: 0.380223 1.20944 2.81665 -0.914152 -0.897172 -12.8678
## 4: 115425.25: 0.380202 1.20948 2.81664 -0.914016 -0.897169 -12.8677
## 5: 115425.25: 0.380180 1.20942 2.81664 -0.913895 -0.897180 -12.8677
## 6: 115425.25: 0.380213 1.20944 2.81663 -0.913896 -0.897239 -12.8675
## 7: 115425.25: 0.380182 1.20946 2.81661 -0.913863 -0.897752 -12.8662
## 8: 115425.25: 0.380213 1.20936 2.81660 -0.913920 -0.898306 -12.8648
## 9: 115425.25: 0.380199 1.20934 2.81660 -0.913717 -0.898888 -12.8635
##
## **Iteration 2
## LME step: Loglik: -20594.53, nlminb iterations: 9
## reStruct parameters:
## subject_id1 subject_id2 subject_id3 subject_id4 subject_id5 subject_id6
## 0.3801986 1.2093374 2.8165952 -0.9137169 -0.8988876 -12.8634672
## Beginning PNLS step: .. completed fit_nlme() step.
## PNLS step: RSS = 5412.717
## fixed effects: 0.8425702 0.9113531 0.9825489
## iterations: 1
## Convergence crit. (must all become <= tolerance = 1e-05):
## fixed reStruct
## 0.000000000 0.002561383
## 0: 115425.25: 0.380186 1.20932 2.81660 -0.913585 -0.899474 -12.8616
## 1: 115425.25: 0.380186 1.20932 2.81660 -0.913585 -0.899474 -12.8616
## 2: 115425.25: 0.380186 1.20932 2.81660 -0.913585 -0.899474 -12.8616
##
## **Iteration 3
## LME step: Loglik: -20594.53, nlminb iterations: 2
## reStruct parameters:
## subject_id1 subject_id2 subject_id3 subject_id4 subject_id5 subject_id6
## 0.3801860 1.2093186 2.8165977 -0.9135847 -0.8994736 -12.8615728
## Beginning PNLS step: .. completed fit_nlme() step.
## PNLS step: RSS = 5412.717
## fixed effects: 0.8425702 0.9113531 0.9825489
## iterations: 1
## Convergence crit. (must all become <= tolerance = 1e-05):
## fixed reStruct
## 0.000000e+00 1.495554e-07summary(mctb_nlme.nlslist_a)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: sooner_choice_scaled ~ (endowment_later_scaled * (b^t0 * d^k_w * pratio)^(1/(a - 1)))/(1 + pratio * (b^t0 * d^k_w * pratio)^(1/(a - 1)))
## Data: select(mctb_data, t0, k_w, pratio, subject_id, sooner_choice_scaled, endowment_later_scaled)
## AIC BIC logLik
## 41209.05 41290.67 -20594.53
##
## Random effects:
## Formula: list(a ~ 1, b ~ 1, d ~ 1)
## Level: subject_id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## a 0.31186232 a b
## b 0.16056237 0.352
## d 0.02550802 0.246 0.610
## Residual 0.42649076
##
## Fixed effects: list(a ~ 1, b ~ 1, d ~ 1)
## Value Std.Error DF t-value p-value
## a 0.8425702 0.009861783 24815 85.4379 0
## b 0.9113531 0.005396787 24815 168.8696 0
## d 0.9825489 0.000824174 24815 1192.1624 0
## Correlation:
## a b
## b 0.331
## d 0.248 0.541
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.644619e+00 -1.424034e-01 -7.529998e-20 4.850362e-01 5.602953e+00
##
## Number of Observations: 25896
## Number of Groups: 1079summary(mctb_nlme.nlslist)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: sooner_choice_scaled ~ (endowment_later_scaled * (b^t0 * d^k_w * pratio)^(1/(a - 1)))/(1 + pratio * (b^t0 * d^k_w * pratio)^(1/(a - 1)))
## Data: select(mctb_data, t0, k_w, pratio, subject_id, sooner_choice_scaled, endowment_later_scaled)
## AIC BIC logLik
## 51574.3 51631.44 -25780.15
##
## Random effects:
## Formula: list(b ~ 1, d ~ 1)
## Level: subject_id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## b 0.13088065 b
## d 0.02301261 0.704
## Residual 0.57988744
##
## Fixed effects: list(a ~ 1, b ~ 1, d ~ 1)
## Value Std.Error DF t-value p-value
## a 0.9454517 0.000794510 24815 1189.9811 0
## b 0.9190270 0.004295194 24815 213.9663 0
## d 0.9829883 0.000715799 24815 1373.2737 0
## Correlation:
## a b
## b -0.027
## d -0.016 0.621
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.578681e+00 -2.167673e-01 -1.581553e-05 2.787358e-01 4.000113e+00
##
## Number of Observations: 25896
## Number of Groups: 1079plot(coef(mctb_nlme.nlslist_a)[,'b'], coef(mctb_nlme.nlslist)[,'b'])
plot(coef(mctb_nlme.nlslist_a)[,'d'], coef(mctb_nlme.nlslist)[,'d'])